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Abstract. 

We investigate the particle and kinetic-energy densities for a system of N fermions 
confined in a potential V^(r). In an earlier paper [J. Phys. A: Math. Gen. 36, 1111 
(2003)], some exact and asymptotic relations involving the particle density and the 
kinetic-energy density locally, i.e. at any given point r, were derived for isotropic 
harmonic oscillators in arbitrary dimensions. In this paper we show that these local 
virial theorems (LVT) also hold exactly for linear potentials in arbitrary dimensions and 
for the one-dimensional box. We also investigate the validity of these LVTs when they 
are applied to arbitrary smooth potentials. We formulate generalized LVTs that are 
suggested by a scmiclassical theory which relates the density oscillations to the closed 
non-periodic orbits of the classical system. We test the validity of these generalized 
theorems numerically for various local potentials. Although they formally are only 
valid asymptotically for large particle numbers iV, we show that they practically are 
surprisingly accurate also for moderate values of N. 
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1. Introduction 

The virial theorem for a particle bound in a local potential ^(r) relates its kinetic and 
potential energies through the following general relation (which we may quote without 
need of referring to any of the standard text books): 

{T) = {r.VV). (1) 

Classically, the brackets (. . .) imply an average over the space covered by the particle. 
Quantum-mechanically, they indicate the expectation values of the corresponding 
operators in a given (eigen-)state of the particle. For a spherical potential homogeneous 
in r, the r.h.s. of ([1]) is proportional to the average potential energy (V); for any other 
different iable V{r) the result is not proportional to, but still an energy related to the 
particle's potential energy, while ( T ) always is the average kinetic energy. 

An essential aspect of the virial theorem ([T]) is that it relates integrated energies to 
each other, averaged over all possible locations of the particle. In the present paper, we 
address the question to which extent a relation (or relations) may be established between 
the kinetic and potential energies locally at any given point r in space. Quantum- 
mechanically, we shall study relations between the corresponding spatial densities, i.e., 
the particle, potential-energy and kinetic-energy densities, valid at any point r. Such 
relations shall be termed here local virial theorems. The systems we are investigating 
consist of fermions bound in a local potential V^(r), and we shall study relations 
between their exact (quantum-mechanical) spatial densities. Although we treat the 
particles as non-interacting, we keep in mind that a local potential V{r) may well 
represent the self-consistent ('mean-field') potential of an interacting system in the 
mean-field approximation, as obtained in the framework of density functional theory 
(DFT) (see, e.g., P). 

Recent experimental success confining fermion gases in magnetic traps [2] has led 
to a renewed interest in theoretical studies of confined degenerate fermion systems at 
zero PSlElEllTlElElITnillllllS] and finite temperatures [El [H]. Quite some effort 
has been devoted in these articles to establish local virial theorems for various types 
of confining potentials. In [HI [H], exact local virial theorems have been established 
for fermions bound in isotropic harmonic oscillator (IHO) potentials in arbitrary space 
dimension D. Some alternative virial theorems involving differentiation or integration 
of their particle density were also given in [H]. Our aim here is to investigate to what 
extent the results of [TT|[n] may be generalised to arbitrary local potentials V{r). While 
an obvious attempt is to simply replace the IHO potential V{r) = cr^ by an arbitrarily 
chosen local potential V^(r) in all those relations, we can only show that this leads to 
exact results for the D-dimensional linear potential V{r) = a. - r with a constant vector 
a (which is not confining, but whose densities can nevertheless be calculated). For 
other potentials we find, however, that the local virial theorems and other relations are 
fulfilled approximately in the limit of large particle numbers N. Formal support of this 
finding comes from a semi classical theory developed recently [151 [HI [13 , in which the 
oscillating parts of the spatial densities are expressed in terms of the closed orbits of 
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the classical system. From this approach, one finds immediately a differential form of 
the basic local virial theorem, stated in Eq. f lTTj) below, which is valid for arbitrary local 
potentials. Our present investigations will therefore be guided to an important degree 
by the semiclassical theory and the understanding of the density oscillations emerging 
from it. 

Our paper is organised as follows. In section |2] we give the basic definitions of 
the quantum-mechanical spatial densities. In section [3] we present analytical results, 
both exact quantum- mechanical ones and their asymptotic limits for — )■ oo, for some 
specific systems: (1) (IHO) potentials and (2) linear potentials, both for arbitrary D 
dimensions, and (3) the one- dimensional box (or infinite square- well potential). We also 
give the Thomas-Fermi (TF) results for the asymptotic average parts of the densities 
and characterize two types of density oscillations that occur for all potentials with 
spherical symmetry in D > 1 dimensions except for IHO potentials. Our generalised 
local virial theorems are then formulated in section HJ after sketching the semiclassical 
theory guiding us to them, and tested numerically for spherical and non-spherical 
quartic potentials and for the two-dimensional circular billiard. Section [5] contains a 
summary and conclusions. Some detailed formulae for linear potentials and for the 
one-dimensional box are given in appendices A and B, respectively, and some (integro-) 
differential equations for the density are briefly discussed in appendix C. 

2. Basic quantum- mechanical definitions 

Let us recall some basic quantum-mechanical deflnitions, using the same notation as 
in [11]. We start from the stationary Schrodinger equation for particles with mass m, 
bound by a local potential V^(r) with a discrete energy spectrum {En}'- 



The potential V{r) can be considered to represent the self-consistent mean fleld of an 
interacting system of fermions obtained in the DFT approach. The single-particle wave 
functions 0n(i") are then the Kohn-Sham orbitals [18] and p(r) is (ideally) the ground- 
state particle density of the interacting system [T9] . 

We order the spectrum and choose the energy scale such that < Ei < E2 < . . . < 
En < ■ ■ ■■ We consider a system with an even number of fermions with spin s = 1/2 
fllling the lowest levels, and define the particle density by 



Here A is the Fermi energy and the factor 2 accounts for the fact that due to spin 
degeneracy, each state is at least two-fold degenerate. Further degeneracies, which may 
arise for D > 1, will not be spelled out but included in the summations over n. For the 




(2) 




n(E„<X) 



(3) 
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kinetic-energy density, we consider two different but equivalent definitions 

' - -^^H C(r)VV™(r) , (4) 

n(E„<A) 

- ^2Y: |V0„(r)r, (5) 

n{En<X) 

which upon integration yield the exact total kinetic energy. Due to the assumed time- 
reversal symmetry, the above two functions are related by 

r(r) = n(r)-i^VV(r). (6) 
An interesting, and for the following discussion convenient quantity is their average 

:= ^[r(r) + n(r)]. (7) 
We can express r(r) and ri(r) in terms of ^(r) and V^p(r): 

r{r) =^W-^|^VV(r), (8) 

n(r)=e(r) + ^^VV(r), (9) 

so that p(r) and ^(r) can be considered as the basic particle and kinetic densities 
characterising our systems. Eqs. ([3]) ^ © are exact for arbitrary potentials V{r). For 
any even number of particles they can be computed once the quantum-mechanical 
wave functions </'ri(r) are known. 

For harmonic oscillators it has been observed long ago [211 [22] that inside the 
system (i.e., sufficiently far from the surface region), ^(r) is a smooth function of the 
coordinates, whereas r(r) and Ti(r), like the density p(r), exhibit characteristic shell 
oscillations that are opposite in phase for r and ti. It has also been noted that ^(r) 
results from the momentum space average of the classical kinetic energy over the Wigner 
transform of the density matrix [221 123] . 



3. Exact and asymptotic quantum-mechanical results 

In this section we first recall in l3. li the results of the Thomas- Fermi theory for the smooth 
parts of the spatial densities which hold for arbitrary local potentials. We then discuss 
exact quantum-mechanical expressions and relations amongst the densities, and their 
asymptotic forms in the limit of large particle numbers N, for some specific potentials. 
In section 13.21 we review known results [TTl [H] for isotropic harmonic oscillators with 
filled shells in arbitrary space dimensions D. In section 13.31 we present new results for 
linear potentials and in section 13.41 for the one- dimensional box with infinitely steep 
walls. In section 13. 5[ finally, we discuss the separation of the spatial densities into 
smooth and oscillating parts and point out the existence of two kinds of oscillations for 
potentials in D > 1 dimensions with spherical symmetry. 
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3.1. Thomas-Fermi limits 

In the limit N ^ oo, the spatial densities are expected to go over into the approximations 
obtained in the Thomas- Fermi (TF) theory [21] . These are given, for any local potential 
V{r), by 

' ' fiflAx.-i'W)"'-'. (11) 



PTF(rj 
TTpfr) 



iD + 2) T{D/2) \2'Kh? 



and 



eTF(r) = (ri)TF(r) = rTF(r) . (12) 

The Fermi energy Atf is defined such as to yield the correct particle number upon 
integration of ptf(i') over all space. The TF densities are valid only in the classically 
allowed regions limited by the classical turning points r^, defined by V{y\) = Atfj so 
that Atf > V{r). Outside these regions the TF densities must be put equal to zero. 
The direct proof that the quantum-mechanical densities, as defined in section [2] in terms 
of the wave functions, reach their above TF limits for — )■ oo is by no means trivial. 
For IHOs it has been given in [TT]. For other potentials, it follows implicitly from our 
results in section HI 

The TF densities (|TO|) - ([TT]) fulfill the following functional relation: 



AttD 

Ttf{v) = rTF[pTF(r)J 



'^''pi;^/^(r), (13) 



2m(D + 2)[4 V2 

which in [17] has been shown to hold also between the exact densities p(r) and r(r) to 
leading order in their oscillating parts. 

For smooth potentials in D > 1 dimensions, next-to-leading order terms in 1/A^ 
modify the smooth parts of the spatial densities. These are obtained in the extended 
Thomas-Fermi (ETF) model as corrections of higher order in h through an expansion 
in terms of gradients of the potential [2S]. These corrections usually diverge at the 
classical turning points and can only be used in the interior of the system, sufficiently 
far away from the turning points. We do not reproduce the explicit expressions of the 
ETF densities here, but refer to chapter 4 of Ref. [26] where they are given for arbitrary 
smooth potentials in D = 2 and 3 dimensions, and to [TT] where explicit results are 
given for spherical harmonic oscillators in D = 2 and 4 dimensions. 



3.2. Isotropic harmonic oscillator in D dimensions 

We review here some exact expressions [HI [H] for the densities in the isotropic harmonic 
oscillator (IHO) potential in D dimensions defined as 

771/ 

V{r) = —u^r^, r=\r\, r G M^, (14) 
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and some equations relating them [TTl [T3], which serve as starting points for our later 
investigations. The eigenenergies En and their degeneracies dn are given by 

En = hj{n + D/2), dn=i ], (15) 

where n = 0, 1, . . . is the principle quantum number. We choose the particle number 
such that the first M + 1 degenerate shells are completely filled, where M is the principle 
quantum number of the last occupied shell, and the densities become spherical. The 
number of particles then becomes 

-(-)-^- (-) 

From some simple expressions for the densities p(r) and ^(r) given in |5], the 
following relation has been shown in [H] to be exact for IHOs with M filled shells: 

where V{r) is given in ( IT4|) . Here A denotes the radial part of the Laplacian operator 
in D dimensions 

^ _ d2 {D-l) d 
dr^ r dr ' 



and Am is defined as 

Am = ^ 



M+-iD 
2 ^ 



(19) 



which corresponds to the mean of the highest occupied and the lowest unoccupied level 
and can be identified with the Fermi energy at zero temperature. 

Since (|T7I) relates the kinetic-energy density ^(r) with the potential-energy density 
V{r)p{r), it represents one form of a local virial theorem, although it involves a term 
proportional to the Laplacian of the particle density. We may eliminate this term in 
favour of the kinetic-energy density T(r), using the relation ([8]), to obtain the relation 

T{r) = [\M-V{r)]p{T)-^i{r). (20) 

In the following, this relation shall be called the basic local virial theorem (LVT), and 
its validity for other than IHO potentials will be investigated. 

Another type of virial theorem, which involves an integral over the density p(r) 
over the whole space, was derived in [1] 



e(r) = - V'{q)p{q)dq, (21) 

where V'{r) = dV{r)/dr is the radial derivative of the IHO potential ( [T4|) . We will in 
the following call (12T]) the semi-local virial theorem (SLVT), since it holds locally for 
the kinetic-energy density ^(r) but requires the knowledge of the density p{r) over the 
whole space. 

All of the above equations are so far known to be exact only if V{r) is the IHO 
potential ([Hj) with M + 1 filled degenerate shells, and if Am is given by ( |T9l) . Their 
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forms, however, suggest immediate generalisations to arbitrary potentials V{r). This is 
one of the main goals of the present paper. 

Other interesting aspects are related to the quantum shell oscillations in the 
densities p(r) and ^(r) which were decomposed into smooth and oscillating terms in 
pT| (see there for the precise definition of the smooth terms) by writing 

p(r) = p(r) + 6p{r) , ^(r) = ^(r) + 6^{r) , 

T{r) = T{r) + 6T{r) , ri(r) = ri(r) + 5ri(r) . (22) 

The following asymptotic behaviours of these quantities were derived in [TT] from an 
expansion of the exact densities in powers of M~^. 

a) In the limit — ?■ oo, the smooth parts of the densities go over into their Thomas- 
Fermi (TF) expressions ffTOj) - f|T2|) given in section ISTT] (or their extensions for D > 1), 
except in a narrow region close to the classical turning points. In the same limit, one 
finds Am — )■ Atf- 

b) The oscillating parts Sp{r), (5r(r) and STi{r) are of order relative to their smooth 
parts, while 6^{r) is of relative order M~^. Practically, 6^{r) can be neglected in the 
interior of the system and ^{r) is essentially smooth there, as observed numerically 
frn. [2Tj . Only close to the classical turning point, S$,{r) becomes comparable in 
amplitude to Sp{r) and 5r(r). 

c) As a consequence of the fact that ^(r) is smooth in the interior of the system, the 
asymptotically leading oscillations in the two kinetic-energy densities r(r) and Ti(r) are, 
due to ([8]) and equal in magnitude but opposite in phase: 

5ras(r) = -(5ri)as(r). (23) 

Here the subscript 'as' refers to the asymptotic large- (or large-M) limit. Deviations 
from this asymptotic relation occur only near the classical turning points. 

d) Extracting from ( 120|) the oscillating terms and neglecting (5^(r), one obtains the 
asymptotic relation 

SrUr) ^ [XM-V{r)]SpUr), (24) 

which we will call the basic differential LVT for the asymptotically leading oscillating 
terms in r(r) and p(r). In fact, this is the form of the LVT that could be derived 
from the semiclassical theory in [151 IE] for arbitrary (also non-spherical) potentials [cf. 
equation ( !7T|) in section 14.11 below] . 

e) For not too large distances r from the centre, the oscillating part Sp{r) is 
asymptotically (up to terms of order M~^) given by 

^-«-(-)"(i^)(B^)"^"'^)' 

where Ju{z) are the standard Bessel functions, and the dimensionless quantities u and 
z are defined by 

u = D/2-l, z = 2rpx/h, px = ^/^^^i. (26) 
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Px being the classical Fermi momentum. The function in f l25p is actually an eigenfunction 
of the kinetic energy operator with eigenvalue 4Am: 

-^A5pUr) = ^XMSpUr). (27) 
2m 

In [11] it was also shown analytically that the asymptotic relation (valid for M — oo) 
6TUr)^XM6pUr) (28) 
is well fulfilled in the interior of the system where the potential can be neglected. 
However, the full differential LVT including the potential holds equally well also 
at larger distances except close to the turning points. This is shown in figure [1] for 




2 4 J. 6 8 10 



Figure 1. Test of the asymptotie relation (HH) for N = 79 422 particles (M = 60) in 
the 3D IHO. Solid line: l.h.s., dotted line: r.h.s. of (p4)) . (Units: h = oj = m = 1.) 

N = 79 422 particles (corresponding to M = 60) in a harmonic oscillator in D = 3 
dimensions. The agreement between the two sides is clearly superior to that obtained 
for f l28|l in |Tl] (see Fig. 5 there). The divergence at the classical turning point is due 
to the ETF correction included in the smooth density T{r). Note that for small r, 
the oscillations are accurately described by ( 125|) which for D = 3 (z/ = 1/2) becomes 
proportional to the spherical Bessel function jo{z). 

f) The TF functional relation (fT3|) was analytically shown (in the limit M oo) to be 
valid also between the exact densities r(r) and p(r) to leading order in 1/M: 

r(r) = rTF[p(r)] + 0(M-2), (29) 

i.e., including the terms 5ras(r) and Spi^{r) which are of order M~^. 

3.3. Linear potential in D dimensions 

After reviewing earlier results for IHO potentials in the previous section, we now present 
new results for a linear potential in D dimensions 

y(r)=a-r, (30) 

with a vector of D constants which we, without loss of generality, assume to be 
positive: 

a = (ai, a2, . . . , od) , ai > . (31) 
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This potential does not bind, but it confines a particle to the left half of space bounded 
by a fiat hyper-surface, leading to a continuous quantum energy spectrum. However, 
this system is of interest, because it allows us to study density oscillations in the vicinity 
of a (more or less steep) surface, the so-called Friedel oscillations (see also appendix A), 
and to regularise a divergence problem in the semiclassical theory (see [T6l fT7] ). 

Since the potential f l30|) is separable, the Schrodinger equation reduces to the one- 
dimensional case of the linear ramp whose solutions are given in terms of Airy functions, 
as is well known from WKB theory [27]. To derive the spatial densities, we start from 
the non-diagonal Bloch density which quantum-mechanically is defined in terms of the 
solutions of ([2]) by 



C(r,r';/3) = J]0:(rO 



r e 



(32) 



where the sum is over the complete spectrum and /3 is a complex variable. Using centre- 
of-mass and relative coordinates q = (r + r')/2 and s = r — r', we may express the Bloch 
density as a function of the variables q, s and /3. The densities p(r) and ^(r) are given 
by the following inverse Laplace transforms of C(q, s; /3) (see, e.g., [26]): 

^ {C(q,s;/3)}q=r,s=ol , (33) 



p(r) 



P 



and 



-{V^C(q,s;/3)}q=r,s=o 



For the linear potential ( l30l) the Bloch density is exactly known (see, e.g., [28] 

\D/2 



C(q,s;/3) 



271^2/3 ) 



(34) 



(35) 



where 



|ap. The particle density then becomes [28] a convolution integral 

(36) 



p(r) = 22/V f pTF(r; X - E) ki{-2'^/^aE) dE . 

J —oo 



where Ptf('"; Atf) is the TP density given in ffTOl) evaluated in terms of the potential 
(15U]) . Ai{z) is the Airy function [21] and a is given by 

1/3 



2m 



(37) 



Performing the derivatives occurring in ( 134|) with the explicit form of (135|) . using 
(Vs- s) = we find 



D 



^C(q,s;/3)}q=r,s=o 



D 
~2 



p(r,A')dA', 



(38) 



whereby the second step is due to a known property of the Laplace transform [29] given 
in ( 136|) . Alternatively, this density can also be written as a convolution integral 

(39) 



^(r) = 22/V I rTF(r; \-E) K\{-2^l'^oE) dE . 

J —oo 
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The proof is easily found by differentiating equations and with respect to A and 
noting from and that drTF(r; A)/dA = (D/2)pTF(r; A) and rTF(r; A = 0) = 0. 

In the results above the Fermi energy A is a continuous parameter, reflecting the fact 
that the spectrum of the potential (l30!l forms a continuum. For this reason, the densities 
( l36l) and (l39ll cannot be normalised. They diverge, in fact, to the far left of the turning 
point. However, we can extract their oscillating parts which will be significant in the 
vicinity of the turning point. As shown in the appendix A, the asymptotic expansion 
of the Airy functions allows us to separate the densities as in (l22i) into smooth and 
oscillating parts. The smooth parts are found to be exactly the TF densities given in 
( ITOl) - ( fT2l) for D = 1, and their ETF extensions [26] for D > 1, while the oscillating 
parts are explicitly given in the appendix A. 

The integrals in ( !36|) and ( l38l) cannot be easily done for arbitrary D. Without 
knowing their explicit forms we can, however, derive the relation (1171) . reading here 

e(r) = 777^ I ^ V V(r) + p(r) [A - ^ (r)] 1 , (40) 



-VV(r) = 22/V [ (-E)pTF(r; \ - E) Ai(-22/VE) dE. (42) 



{D + 2) {8m 

whereby V^(r) now is given by (1301) . To prove it, we first use the identity 

VVTF(r, X-E) = pTF(r, A - E) , (41) 

which holds for the potential (130|) . under the integral of ( 136|) . perform two integrations 
by parts and use the differential equation [29] Ai"{z) = zAi{z) and (137|) to find 

8m 

Combining now the three terms in the square brackets on the r.h.s. of ( 140|) before 
integrating and using (15^ and fl42p . the integrand becomes, apart from the factor 
Ai(-22/VE) 

[A - V{r) - E] pTF(r; \ - E) = ^^^^ rTF(r; \ - E) , (43) 

which with (139|) leads directly to (140|) . Using the same manipulations as in section I3l2| 
we find the LVT given in (!20|) . 

In the appendix A we show that the SLVT fl2T|) is valid also for the linear potential 
fl30|) exactly for D = 1, as given in flA.6p . For arbitrary D > 1, one may formally 
write the density as a multiple convolution integral of D one-dimensional densities 
of the form (lA.Sp . because the D-dimensional Bloch density (!35|) is a product of D 
one-dimensional Bloch densities. Unfortunately, these convolution integrals can not 
be done analytically. However, explicit results can be found if one restricts oneself to 
projections of the densities along an arbitrary Cartesian axis Xi {1 < i < D), so that 
r = (0, . . . , Xj, . . . , 0). For D = 1 this is, of course an exact result. For the present, we 
use the simplified notation 

p(a;,) = p(0,...,0,a;„0,...,0), (44) 
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and likewise for the other densities. Along the Xj axis, the density ( I36p is only a function 
of aiXi — A, so that the integral in f p8|l can be performed as in the one-dimensional case, 
yielding the generalisation of ( 1A.6I) : 

D f°° 

= j J aip{x[) dx[ . (45) 

This expression is identical with the SLVT ( I2T]) for the IHO potential in D dimensions, 
when the radial variable r there is replaced by the coordinate Xi and the potential fl30l) 
is used. 

We have thus found the interesting result that for the linear potential ( 130|) in D 
dimensions, the spatial densities along any Cartesian axis fulfill the same local virial 
theorems as for the IHO potentials in D dimensions. Note that for the IHOs they only 
hold for the specific values f|T9|) of A^/. In the present case, however, they are valid for 
arbitrary values of A, since there is no shell structure in the continuous energy spectrum 
{E} of the linear potential ( 130|) and A is a smooth function of the energy E. 



3.4- The one- dimensional box 

Another system, for which the wave functions are known analytically, is the one- 
dimensional box with length L and ideally reflecting walls (corresponding to Dirichlet 
boundary conditions for the wave functions, see appendix B): 

V{x) = for < X < L , V{x) = oo else . (46) 

Detailed calculations for the densities are given in appendix B. It suffices here to state 
the main results regarding the local virial theorems. The oscillating part of the density 
asymptotically satisfies the relation 

-^6p'Ux) = AXr^6p{x), (47) 

This is the equivalent of ( 1271) valid asymptotically for IHOs. It is also easy to show 
that the differential LVT (!24|) derived for IHOs is satisfied here, too, with the proviso 
V{x) = inside the box: 

6t^s{x) = \tfSp{x) . (48) 

Furthermore, as shown in appendix B, the oscillating parts of the two forms of kinetic- 
energy density also fulfill the relation 

6ti{x) = -6t{x) . (49) 
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3. 5. Structure of the oscillating parts of the densities in radial potentials 

Based on the results discussed above, the spatial densities may be decomposed in the 
following way: 

P(r) = P(E)TF(r) + 5p(r), (50) 

r(r) =r(E)TF(r) + 5r(r), (51) 

ri(r) = (ri)(E)TF(r) + <5ri(r), (52) 

e(r) =e(E)TF(r) + 5e(r). (53) 

For one- dimensional systems and for billiards in arbitrary dimension D, the subscripts 
TF hold and hence the explicit relations ( fTOj) - ( fTT]) can be used [30]. The oscillating 
parts, denoted by the symbol 6, have been approximated semiclassically in [151 [HI [H] 
as discussed in section 14.11 below. 

The systems discussed above in this section are the only ones, to our knowledge, 
in which explicit expressions for the oscillating parts of the spatial densities can be 
extracted. Numerically, however, we have investigated the densities in several potentials 
in D > 1 dimensions with radial symmetry such that V{r) = V{r), where r = |r|. 
We have observed that the function ,^(r) for D > 1 in general is not smooth in the 
interior and does not therefore coincide asymptotically with the corresponding (E)TF 
approximation, such as is the case for isotropic harmonic oscillators. Indeed we find 
that ^(r) contains oscillations whose amplitudes are comparable to - and in higher 
dimensions D > 2 even larger than - those of the regular fast shell oscillations appearing 
in the densities p(r), r(r) and Ti(r) for harmonic oscillators. They are, however, rather 
irregular and have a longer wave length in the radial variable r. 

An example is shown in figure |2] for a spherical billiard with unit radius containing 
= 100068 particles. Note the irregular, long-ranged oscillations of ^(r) around its 
bulk value [30] Ctf seen in the upper panel. In the lower panel, where we exhibit only 
an enlarged region around the bulk value, we see that r(r) and Ti(r) oscillate regularly 
around ^{r), but much faster than ,^(r) itself and with opposite phases. The same two 
types of oscillations are also found in the particle density p(r). 

For radial systems, we can thus decompose the oscillating parts of the spatial 
densities defined in (!50|) - ( !53|) as follows: 

6p{r) = 5rP(r) + 5i„pW , (54) 

5r(r) = 6rT{r) + 6i„T{r) , (55) 

(5ri(r) = 6rTi{r) + 6i„Ti{r) , (56) 

=SiMr). (57) 

Here the subscript "r" denotes the regular, short-ranged parts of the oscillations, while 
their long-ranged, irregular parts are denoted by the subscript "irr" . We emphasise that 
this separation of the oscillating parts does not hold close to the classical turning points. 

As we see in figure [2] and in later examples, the oscillating parts defined above 
fulfill the following properties in the interior of the system (i.e., except for a small 
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Figure 2. Kinetic-energy density profiles of a 3-D spherical billiard with N = 100068 
particles (units: lT?/2m = i? = 1). Upper panel: ^(r) (solid line) and its constant TF 
value ^TF (dashed). Lower panel: T{r) (dashed), Ti(r) (dotted) and ^(r) (sohd line). 
Note that in both panels, the vertical scale does not start at zero. 



region around the classical turning points): 

a) For D > 1, the irregular oscillating parts of r(r) and ri(r) are asymptotically identical 
and equal to S^{r): 

5i„r(r) ^ 5i„ri(r) ^ S.^r) = ^r) ■ (58) 

b) The irregular oscillations are absent (i.e., asymptotically zero) in the densities of all 
potentials in D = 1 and, in addition, in the IHOs ( |T4l) and the linear potential (130|) for 
arbitrary D. 

c) The regular oscillating parts of r(r) and ri(r) are asymptotically equal with opposite 
sign: 

S,T{r) ^ -6,Ti{r). (59) 

This relation holds in particular for harmonic oscillators for which it has been derived 
in [11], as given in (123|) . 

All these properties could be explained by the semiclassical theory developed in 
[T5l [T7] , whose main results will be summarized in section 14.11 below. We anticipate 
here that the fast regular oscillations are due to linear radial (i.e., self-retracing) classical 
orbits, while the irregular slow oscillations are due to non-radial (i.e., not self-retracing) 
classical orbits. The oscillations in ^(r), however, are due only to non-radial orbits (if 
they exist) and are therefore of the irregular type. 

In the following, the symbol S denotes the sum of both types of oscillating parts; 
the subscripts will only be used if reference is made to one particular type of oscillations. 
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4. Generalised local virial theorems 

So far we have presented exact local virial theorems (LVTs) that were derived purely 
quantum mechanically. They were shown in section [3] to hold both for IHOs [11] and for 
linear potentials in arbitrary dimension D (for the latter along any of the Cartesian 
coordinates). Some asymptotic relations for the oscillating parts of the quantum- 
mechanical densities have been given, too, and shown to hold also in the one-dimensional 
infinite square well. 

In the present section we shall investigate to what extent these relations can be 
generalized to arbitrary different iable local potentials V^(r). Since we have no exact 
proofs except for the potentials mentioned above, we employ a semiclassical theory of 
density oscillations developed recently in [151 [El IH] • This theory is asymptotically valid 
in the limit ^ — which, for the systems under investigation here, corresponds to the 
limit — 7- oo. The equivalence of these two limits can directly be seen from equation 
( 1T9|) for spherical harmonic oscillators. For arbitrary local potentials, it follows from 
the general validity of semiclassical quantization in the limit of large quantum numbers 
which, for finite classical actions, is the same as the limit /I — )■ (see, e.g., Ref. |26]). 

Correspondingly, the generalized virial theorems presented below are not exact, but 
asymptotic theorems that are expected to apply for large particle numbers A^. As we 
will see, however, they work also surprisingly well for moderate values of A^. 

We will briefiy sketch the semiclassical theory in section |4?T1 and in sections 14.21 and 
14.31 we shall present the generalized LVTs and test them numerically for some specific 
potentials. 

4.I. Sketch of semiclassical theory for density oscillations 

We reproduce here the main formulae for the semiclassical approximations to the 
oscillating parts of the spatial densities, which were derived in [151 IIHl [IZ] from 
the semiclassical Green function established by Gutzwiller [311 132] . Starting from 
the decompositions f l50|) -f l53|) . the following expression for the oscillating parts of the 
densities are valid to leading order in h: 

6p{r) ^ Re a, A,CK r) e*^^^''') , (60) 

7rp(A,r) ^ 

Sr{r) ^ Re a„ V ^,(A, r) e*-(^'^) , (61) 

7 

5ri(r) ^ Re y Q (A, r) ^^(A, r) e^^-^^'"-) . (62) 



TT 

7 



The sums are over all orbits 7 of the classical system that lead from a point r back to 
the same point r. The phase function $^(A,r) is given by 

<l>^(A,r) = 5'^(A,r,r)/;i-/i^^, (63) 
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in terms of the general action integral along the orbit 7, taken at the smooth (ETF) 
Fermi energy A = A(e)tf 



S,{X,r,r') = J p(A,q)-dq, (64) 
where p(A, r) is the classical Fermi momentum 



r 



p(A, r) = ^\/2m[X - V{r)] , p{\, r) = |p(A, r) | , (65) 



defined only inside the classically allowed region where A > V^(r). The Morse index 
/i^ is equal to the number of conjugate points along the orbit [32]. The semiclassical 
amplitudes ^-^(A, r) are given by 

AJ\,r) = ^' ■ (66) 

T,(A,r) 

Hereby is the reduced Van Vleck determinant [311 ES] 

= det(9p^/9r'x) , (67) 

where p_L and are the initial momentum and final coordinate, respectively, transverse 
to the orbit 7. T-y(A,r) = dS'^(A, r, r)/dA is the running time of the orbit 7 [33]. The 
"momentum mismatch function" Q^{X,r) appearing in (|62l) is defined as 

g^(A,r) = cos[^(p,p')], (68) 

where p and p' are the short notations for the initial and final momentum, respectively, 
of a given closed orbit 7 at the point r, which are obtained from the action integral ( 164]) 
by the canonical relations 

VrS,(X,r,r') =-p, Vr'^,(A,r,rO = p' . (69) 

r=r' r=r' 

The overall prefactor ao, which depends explicitly on the dimension D, is given by 

an = 27r(2z7r;i)-(^+i)/l (70) 

In principle, all closed classical orbits contribute to the sums in ( I60l) - (l62l) . However, 
as discussed extensively in [17], it is the non-periodic orbits that are responsible for 
the oscillations in the densities. Periodic orbits need to be included in connection with 
uniform approximations necessary at singular points, where the semiclassical amplitudes 
A-y diverge and have to be regularized. (These singular points are the turning points, 
bifurcation points, or r = in systems with radial symmetry; see [T6l [T7] for details.) 
Note that = +1 for p = p', i.e., for periodic orbits, and = —1 for p = — p', i.e., 
for self-retracing non-periodic orbits, in particular for orbits oscillating along a straight 
line which we will call "(radial) linear orbits" below. In one-dimensional systems, there 
are only linear orbits and it could be strictly shown [17] that only the non-periodic orbits 
contribute to the density oscillations. 

It should be stressed that the above expressions do not hold near the classical 
turning points where the amplitudes A.y diverge. They can be regularized by special 
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techniques for which we refer to [T^ . The following relations which we can derive directly 
from these expressions hold therefore only sufficiently far from the turning point. 

Comparing the prefactors in the expressions ( l60i) and ( 16T|) . and using (l65l) . we find 
directly the relation 

5r(r) ^[X- V{r)] 6p{r) . (71) 

This is exactly the differential LVT ( IMl) that was derived [H] for IHOs with M filled 
main shells in the limit M — )■ oo, with the corresponding Fermi energy A = Am given 
in ( lT9ll . In section [3731 we showed it to be fulfilled also for linear potentials at arbitrary 
Fermi energies A. Semiclassically, however, ( I7T]) is valid for arbitrary local potentials 
and arbitrary (even) particle numbers N, since the sums over the orbits 7 cancel from 
( I7T1) and no assumption about the nature of the local potential V{r) has been made at 
this point. The Fermi energy A hereby is that of the (E)TF theory, i.e., A = A(e)tf- 

As shown in [161 IE]) the linear non-periodic orbits always lead to rapid regular 
oscillations 6^p{r) etc., while the irregular oscillations 5i„p{r) etc. are due to the non- 
linear (i.e., more-dimensional) orbits. This explains the observed fact that no irregular 
orbits are found in one-dimensional systems. They are also absent in IHOs and the 
linear potential, since there exist no non-linear non-periodic orbits in these systems; the 
kinetic-energy density ^(r) is therefore smooth and close to ^(e)tf('^) (except possibly 
near the turning points). Looking at the expressions ( |6T1) . ( 162|) and noting that = —1 
for the linear orbits, as stated above, we see immediately that the relation 

5,T{r) = -5,n{r) , (72) 

obtained asymptotically for IHOs, linear potentials and the one- dimensional box in 
section 131 is semiclassically valid for arbitrary potentials V{x) m. D = \, and for arbitrary 
potentials V{r) with radial symmetry 'm. D > \. For the latter one also finds [17] that, 
to leading order in the rapidly oscillating part of the density fulfills the following 
differential equation 

-^^5p,{r) = [\-V{r)]5p,{r). (73) 
8m 



which is the generalization of fl27j) for arbitrary systems with radial symmetry. Close 
to r = where the potential can be neglected, i.e. where V{r) <^ A, fl75|) becomes the 
universal Laplace equation 

-^V^^p{r)^\5,p{r), (74) 
which is the generalization of ( 12 5 p valid for IHOs, with the universal solution 

5,p(r) = (-1)"'—^ (^XU2rp,/h) . (75) 

Here Ju{z) is a Bessel function with index u = D/2 — 1, M -|- 1 is the number of filled 
main shells [M], ^ri is the period of one full radial oscillation and px = (2mA)^/^ is the 
Fermi momentum at r = 0. The expression (175|) was found, indeed, to describe the 
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rapid oscillations of the particle density in spherical potentials (and for D = 1) close to 
the centre very well [151 113 [IS] • 

After compiling these general results derived from the semiclassical approximations 
( 160|) - ( 162|) to the density oscillations, we are now ready to propose the generahzed 
local virial theorems. As already mentioned, the explicit semiclassical expressions given 
above do not apply in the surface regions near the classical turning points without 
additional regularizations [161 HZ]- We therefore will state the theorems below in such 
a way that additional terms, which should only be used in the surface region, appear 
in curly brackets {...}; we shall call them the "surface corrections". Omitting them 
yields the theorems expected to be approximately valid in the interior of the systems. 
Adding them will improve the relations near the classical turning points but may spoil 
their validity in the interior. Furthermore, these surface corrections are only expected 
to be valid for smooth potentials, since they are justified by the local linearization of 
the potential at the classical turning points. 

A rough estimate of the size of the surface region, where these corrections are 
needed, is given by the break-down of the semiclassical approximation near the classical 
turning points. This occurs when the action of the leading closed orbit 7 becomes smaller 
than h, i.e., when S'^(A, r, r) < h. Its precise value depends, of course, on the potential. 
Practically, it is of the order of the wave length h/2px of the Friedel oscillations (see 
appendix A and Ref. [T7]). 

4-2. The local virial theorem (LVT) 

For arbitrary local potentials V{v) in D dimensions, we propose the approximate 
generalized differential LVT: 

5r{v) ^ \\-V{v)]5p{r) [-^5i{v)]. (76) 

The part without the surface correction is just fl7T|) proved semi classically for arbitrary 
potentials. The surface correction is justified by the fact that including it and adding the 
smooth (E)TF densities on both sides leads to the full LVT in ( 120|) . proved for IHOs and 
and shown in section 13.31 to hold also for linear potentials. Since any smooth potential 
can be approximated linearly (or quadratically) near the classical turning points, we 
expect the corrected LVT to be approximately valid in the surface region. 

In order to demonstrate the validity of the differential LVT ( 1761) for a non-spherical 
system, we presently test it for the coupled two-dimensional quartic oscillator 

V{x,y) = \{x' + y')-Kx^y\ (77) 

whose classical dynamics is almost chaotic in the limits n = 1 and k — )■ — cxo [351 136] , 
but in practice also for k = 0.6 (see, e.g., [2Z])- We have computed its wave functions 
using the code developed in [37] . 

In the upper panel of figure [3] we compare left and right sides of (fTUj) for this 
system with N = 632 particles, plotted along the line y = x/y/S as a functions of x. 
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Figure 3. Oscillating parts of spatial densities of = 632 particles in the nearly 
chaotic potential (|77|) with k ~ 0.6 {h ~ m = 1). Top: The solid (black) line gives 
6T{x,y), the dashed (red) line the r.h.s. of the LVT ([76| without surface correction, 
and the dotted (blue) line gives the r.h.s. of (f76| including the surface correction. 
Bottom: SS,{x,y). All results are plotted versus x along the line y = x/\/3. 

The solid line shows the exact 5t{x,x/\/3). The dashed line shows the r.h.s. of the 
LVT ( 176|) without, and the dotted line with the surface correction; both are evaluated 
with the exact 5p{x,x/\/3) and S^{x,x/\/3). We see that the agreement without surface 
correction (dashed line) is very good in the interior; only in the surface region is there a 
visible disagreement. This disagreement is clearly reduced when the surface correction 
is added (dotted line), but at the expense of a less good agreement in the interior. The 
quantity 6^{x, x/y/3) is shown separately in the lower panel and seen not to be negligible 
anywhere. 

Next we generalize the LVT derived in the form ( 120|) for IHOs and shown to be 
valid also for linear potentials. For arbitrary local potentials V{r), we propose the 
approximate generalized LVT: 

r(x) ^ [X-V{x)]p{x)-2^{x), {D = l) (78) 

r(r) ^ [A-V(r)]p(r)--^eETF(r){-|5e(r)}. {D > 1) (79) 

Our justification for this generalization is the following. First we note that the TF 
densities ([TU]) . f ITT]) fulfill exactly the relation 

rTF(r) = [A - V{r)] p^r) - ^ ^TF(r) , (80) 

so that, to leading orders in h, the smooth parts of the relations (ITHl) and fl79l) are exactly 
true. Adding now the differential LVT fl7T]) to the above and using fISUl) - flS51) . we arrive 
at (179ll for D > 1. For D = 1, ^{x) exhibits no oscillations in the interior, so that we 
may add S^{x) everywhere. 
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Figure 4. Test of the generalized LVT ([75)1 = 40 particles in the one-dimensional 
potential V{x) = x^j^. Solid line: exact T(r), crosses: r.h.s. of (j78p using the exact 
densities p(x) and ^(x) (units: h = m ~ 1). 

We first test the one- dimensional LVT (ITHj) in figure H] for the potential V{x) = x^jl 
with = 40 particles. The solid lines show the exact t{x) and the crosses the r.h.s. 
of ( ITHl) . calculated with the exact densities p{x) and ^{x). The agreement is seen to be 
perfect everywhere. 
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Figure 5. Test of the generalized LVT ([7^ without surface correction for N ~ 6956 
particles in the two-dimensional potential V(r) = r^jl. Solid line: exact r(r), dotted 
line: r.h.s. of (|79p using the exact density p{r) (units; h = m = I). 

The LVT (1791) without surface correction is tested similarly in figure [5] for the two- 
dimensional radial potential V{r) = r'^/A with = 16906 particles. The inserts show 
the central and surface regions on enlarged scales. The agreement is again very good; a 
small deviation occurs only near the classical turning point (see the upper right insert) 
where ^etfI?") misses the exponential tail. 
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In figure M we show the same results including the surface correction in f l79p . 
The agreement is now perfect in the surface region; the price paid for this is a slight 
discrepancy near the centre of the system which, however, is not serious. Practically, 
one may therefore live with the surface-corrected LVT ( I7T]) in the whole space. 
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Figure 6. Same as figure [5] but including the surface correction on tlie r.li.s. of ([79|) . 

Although we have in this section restricted ourselves to differentiable potentials, 
we show in figure [7] that the LVT (1791) without surface correction applies also to billiard 
systems. Here we test it for the two-dimensional circular billiard with = 68 particles. 
Close to the surface the LVT does not apply, as expected, but in the interior it works 
surprisingly well even for this relatively small particle number. 
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Figure 7. Test of the LVT fT^l) without surface correction for the two-dimensional 
circular billiard with TV = 68 particles. Solid line: exact T(r), dashed line: Xp{r) with 
A = 160.68303 (units: 11^/2771 = i? = 1). 
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4-3. The semi-local virial theorem (SLVT) 

We next want to generalize the SLVT given in ( l2T|) for IHOs and shown in section 13.31 
to be exact also for linear potentials if r is replaced by any Cartesian coordinate. For 
one-dimensional systems, it can actually be proved to be exact for any differentiable 
potential V{x). It reads 

1 f°° 

^[x) = - V'{x')p{x')dx'. {D = l) (81) 
Taking the derivative on both sides leads to 

ax) = -\v'{x)p{x). (82) 

This equation is easily proved by taking the derivative of the one- dimensional Schro- 
dinger equation ([2]) for each state 0n(a:), multiplying the result by </>Ji(a;) from the left, 
summing over all occupied states up to the Fermi energy and using the definitions of 
the densities. Integrating ( !82l) . noting that the integration constant must be zero since 
all densities vanish exponentially at infinity, leads back to the SLVT flHTl) . (See also the 
discussions in [22l [38| [39].) 

For arbitrary differentiable potentials V{r) m D > \ with radial symmetry, we 
propose the approximate generalized SLVT: 

eETF(r){+5^(r)} ^ ^ j^V'{r')p{r')dr'. {D > I) (83) 

We justify this semiclassically by the following argument. Like above, we note that the 
TF densities flTU]) . f lT^ for spherical potentials fulfill exactly the relation [ID] 

D 

eTF(r) =-J^ V'ir')pTFir')dr'. (84) 

Adding 6p{r') under the integral on the r.h.s. above leads, to leading orders in h, 
to the r.h.s. of (l83l) . However, an integration over the radial variable r' applied to 
the semiclassical expression ( 1601) of 6p{r') yields a factor proportional to h and hence 
suppresses all oscillations in the interior (to leading order in h). This is why only the 
smooth part of ^(r) is contained on the l.h.s. of fl83|) without surface correction. The 
surface correction in ( l83l) leads to the full density ^(r) on the l.h.s. and hence corresponds 
to the SLVT valid exactly for IHOs (!2T|) and linear potentials (H5l) . 

The SLVT without surface correction is tested in figure M for the same system as 
in figure [51 We see that, indeed, the r.h.s. of f[83|l (dotted line) is perfectly smooth and 
can hardly be distinguished from the density ^etf('") (solid line), except very close to 
the surface where the latter lacks the exponential tail. In figure [9] we show the same test 
after adding the surface correction on the l.h.s. of ([83]) . We see that the full ^(r) in the 
interior has the characteristic irregular oscillations which are absent from the integral 
on the r.h.s. of f[83l) . In the surface, however, both sides agree perfectly and have the 
same exponential tail. 
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Figure 8. Test of SLVT ([55]) without surface correction for the same system as in 
figure [5l Solid line: ^etf(^); dotted line: r.h.s. of (|83p 




Figure 9. Same as figure [5] but including the surface correction on the l.h.s. of (|83p . 
The sohd hne here is the full ^(r). 



The integral on the r.h.s. of (1831) is in itself an interesting quantity. Let us call it 
^2{f') by defining, for any dimension D, 

^2(r) = - / V'{r')p{r') dr'. (85) 



2 . , 

Integrating over the whole space in (hyper)spherical coordinates yields 

r) o roo poo 

^2(r)d^r = r^~Ur V'{r')p{r')dr' , (86) 



2 _ 

where fi^ is the integrated solid angle in D dimensions. After integration by parts and 
noting that the densities vanish at infinity, we obtain 

^2(r) d^r = l [ rV'{r)p{r) d^r . (87) 



2 _ 

This is nothing but the r.h.s. of the standard (integrated) virial theorem ([1]) for a 
spherically symmetric potential, and hence identical with the total kinetic energy. Thus, 
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integration of the surface-corrected SLVT (1831) on both sides yields the standard virial 
theorem which is exact. Consequently, the difference between C,{r) and ^2(1^) can be 
written as a local error term R2{r) which integrates to zero and vanishes at infinity: 



As shown in section [3l we know that R2{r) = for IHOs and linear potentials. It 
would be interesting to study mathematically the function R2{r) for other differentiable 
potentials V{r) with spherical symmetry. 

5. Summary and concluding remarks 

This paper deals with local virial theorems (LVTs) that connect kinetic and potential 
energy densities with particle densities for non-interacting fermions, bound in a local 
potential ^(r), at any given point r in space. We have first reviewed exact relations 
that were earlier derived for D-dimensional isotropic harmonic oscillators (IHOs), and 
then proved the same relations to hold also for linear potentials in arbitrary dimensions, 
as well as for the one- dimensional box with Dirichlet boundary conditions. We then 
showed that the LVTs can be generalized to arbitrary local potentials, if they are taken 
as approximate relations, valid asymptotically for large particle numbers N. Practically, 
however, they are found to work numerically quite well also for moderate values of A^. 

Our generalized approximate LVTs are supported by a semiclassical theory, 
developed recently [IHl UHl [IT] and summarized in section 14.11 which relates the 
oscillating parts of the spatial densities to the closed (non-periodic) orbits of the classical 
system. The basic differential LVT ( 1711) was semiclassically shown to hold for arbitrary 
local potentials. It is therefore (asymptotically) valid also for an interacting N-fermion 
system bound by the self-consistent Kohn-Sham potential. We have shown numerically 
that these generalized theorems are well fulfilled for various local potentials. 

Since the semiclassical approximation breaks down at the classical turning points, 
the generalized local virial theorems are not valid in regions close to the surface, roughly 
given by a distance h/2px perpendicular to the closest turning point (where p\ is the 
Fermi momentum). For these regions, we have proposed "surface corrections" to the 
LVTs for smooth potentials that were derived from the local linear approximation to 
the potentials at the turning points and numerically tested successfully. 

We note that, as a direct consequence of the differential LVT ( 17T]) . the TF functional 
relation (fT3|) has been shown in to be valid between the exact densities r(r) and p(r) 
to first order in their oscillating parts for arbitrary local potentials: r(r) ^ rTF[p(i")] 
(except close to the classical turning points). A related result in one dimension, based 
on semiclassical (WKB) arguments, can be found in Ref. [H], where also gradient 
corrections to the TF kinetic energy functional are discussed. 

For systems with spherical symmetry, two kinds of oscillations in the spatial 
densities can be separated, as implied in equations ( 15^ - ( 1571) . In the semiclassical 
theory, the regular, short-ranged ones (denoted by the symbol 6,-) are attributed to linear 
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non-periodic orbits in the radial direction, and the irregular, long-ranged ones (denoted 
by ^ii-r) are due to non-linear orbits and therefore only exist in D > 1 dimensions. This 
also explains the empirical fact that, for all one-dimensional systems and for IHOs and 
linear potentials in any dimension D, the kinetic-energy density ^(r) defined in ([7j) has 
no regular oscillations, since these systems contain no closed non-linear, non-periodic 
orbits. 

An interesting object is the quantity C,2{^) defined in fl85l) . In D = 1 dimension, 
we have shown it to be identical with the exact quantum-mechanical ^(x) for any 
different iable potential V{x). Its identity with ^(r) holds in D > 1 dimensions, too, 
for IHOs and for linear potentials (when taking r to be any Cartesian coordinate), for 
which C,{r) is smooth, as shown in section |3l For arbitrary spherical potentials V{r) in 
D > 1, we expect it to be approximately equal to ^(r) only in the surface region near 
the classical turning points, while in the interior of the systems, it yields the smooth 
part ^etf('^) only, as expressed in the generalized semi- local virial theorem (183|) . 

We expect that our generalized LVTs might be of practical use in the analysis of the 
spatial (kinetic-energy and particle) densities of trapped fermionic atoms. In particular, 
we propose it as a challenge for the cold atoms community to verify the differential LVT 
( I7T]) experimentally. 

In the appendix C we briefly discuss some (integro-) differential equations for the 
particle density p(r) alone, valid in IHOs and linear potentials. Their generalization for 
D > 1 is, however, of little practical use, since it involves also explicitly the regularly 
oscillating part 5rp{r) in the interior of the system, see equation (IC.Sp . which is a priori 
not known. 
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Appendix A. Explicit densities and relations for linear potentials 

In this appendix we give some explicit analytical results for the spatial densities in the 
linear potential (|30|1 in those cases where we have been able to find them. 

Appendix A.l. D = 1 

For D = 1 with V{x) = ax, the expression ( l36i) was found in [28] to be equivalent to 

p{x) = 2^^j' Aiy{ax-E)]dE. (A.l) 

Using the dimensionless variable z defined by 

z = a{ax - E) , (A.2) 
we can rewrite it as 

/•oo 

p{x)=pJ ki\z)dz, (A.3) 

with 

f2maV'^ 

Zx = a{ax - A) , Po = 2 ( j =2aa. (A.4) 



Next, we note [29] that the function w{z) = Ai (z) fulfills the differential equation 
w = w'"/2 — 2zw'. Using this for the integrand of (1A.3P and the differential equation 
for the Airy function as above, we obtain after integration by parts 

p{x) = po{[A^{zx)]'-zxAi\zx)}. (A.5) 
For the kinetic-energy density ^(x) we can rewrite the integral in (138|) for D = 1, using 



as 



1 



oo 



e(x) = -y ap{x')dx'. (A.6) 

This expression is identical with the relation (12T]) obtained for the one-dimensional 
harmonic oscillator {D = 1) when substituting V{x) = ax for the potential. 

As in the case of flA.ip . the integral in (IA.6|) can be done analytically to yield 

e(x) = -^{Ai{zx)Ai'{zx) + 2zx[Ai'{zx)r- 2zlAi\zx)}. (A.7) 

From (]A.5|) we get 



■im 



p"{x) = -aAi(CA)Ai'(CA) , (A.8) 



and using we find 



ix) = y {Ai{zx)A\'{zx) - zx[Ai'{zx)]' + zlAi'izx)}. (A.9) 
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In order to extract the average and leading oscillating components of these densities, 
we use the asymptotic expansion of the Airy function and its derivative [29] for —z ^ 1: 

1 



Ai(-z) 



Ai'(-^) 



^1 11/4 sin(C + 7r/4) - ^ cos(C 



11/4 



TT 



cos(C + 7r/4) 



7ci 
5C 



sin(C 



vr/4) 
-7r/4) 



with 



Cl 



5_ 
72 



C 



1 3/2 



(A.IO) 



(A.ii; 



up to terms of order ( ^. Inserting the above into 
terms up to 0{(~^), we obtain 

p{x) = p{x) + 5pUx) + 0{h) , 

where the smooth part is the TF density 



31) for the density and keeping 

(A.12) 



p{x) 



pTF[X) 



2m 



ax 



(A.13) 



in agreement with ffTOl) . The leading-order oscillating term for ax ^ A simplifies to 

1 1 



5pas(a;) 



cos(2Ca 



(A.14) 



(A.15) 



27r [x — X\) 

with the turning point x\ and the quantity Cx given by 

2 

XA = A/a, = -|z;,|3/2. 

This surprisingly simple-looking expression flA.14p (in view of the complicated nature 
of the Airy function) has a direct semiclassical interpretation in terms of the shortest 
closed classical orbit of the system [T7] . 

Figure lAll shows the exact result (lA.Sp by the solid line. The asymptotic result 
(IA.14P is shown by the dotted line. Although it diverges at the turning point xa, it 
is seen to reproduce the exact 5p{x) even rather close to it. The oscillations, whose 
amplitude reaches a maximum just before the turning point, are the so-called Friedel 
oscillations. 

The oscillating part of ^(x) becomes 
a 1 



sin(2CA) + 0(Ca 



(A.16) 



Note that, since C, oc cr"^/^ oc /i ^, the leading term in 5^{x) is of one order in h higher 
than 5pas{x). Using ([8]), ([9]) and the asymptotic form of flA.Sp 

^p"(x)~f cos(2Ca) + 0(/1), (A.17) 
we find that the oscillating terms of r(x) and ti(x) at the leading order are given by 
5t^{x) = -i5ri)Ux) = cos(2Ca) • (A.18) 

ZTT 
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Figure Al. Oscillating part of spatial density in the one-dimensional linear potential 
pop with a 1, evaluated at the Fermi energy A = 20 (units: h = m = 1). Solid line: 
exact result <\A.5^ : dotted line: asymptotic expression ()A.14|) . 

This is exactly the asymptotic relation ( 123|) obtained for IHOs. Comparing with (]A.14I) . 
we finally get the differential LVT ( 124|) for the linear potential: 

^Tas(a;) = 5pas(a;)(A - ax) , (A. 19) 

valid sufficiently far away from the turning point. 

As discussed in detail in , the closed orbit responsible for the Friedel oscillations 
is the primitive self-retracing orbit (in fTEl [T7] called the "+" orbit) that in general goes 
from a point r to the closest turning point and from there back to r. The wavelength 
of these oscillations near the surface is given by h/2px, where p\ is the classical Fermi 
momentum, as already noted long ago j42j . 

In passing, we note that for the diagonal Bloch density for D = 1, C{x; (3) = C{q = 
x,s = 0] (3) given by ( l35ll . the following differential equation is identically fulfilled: 

|_C(x;/3)-^aC(x;/3)=0. (A.20) 
ox 2 

This is exactly the equivalent of Eq. (A5) given in the appendix of Howard et. al., [10] 
for the harmonic oscillator in D dimensions, but rewritten here for D = 1 and the 
potential ( l30l) (note that the sign in front of the last term in (A5) of [10] is wrong; it 
should be "+"). 

Appendix A. 2. D > 1 along a specific axis 

Specific analytical results can be found for odd values of D. The integral in (136|) for 
D = 3 along the axis Xi can be done by parts, using the explicit forms of the TP density 




5m 



93 



X 



C{x;P)- 



d_ 

dp 



+ ax 
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Pixi) = -7^P%\ Ai{zix)Ai'iza) + 2zix[Ai{zix)f 



I^U^| ioT D = 1 and D = 3, to yield 

1 

-2zlAi'{za)}, {D = 3) (A.21) 

where pio = 2crjaj and the argument zx is given by 

Zix = ai{aiXi - X) , i = l,...,D, (A.22) 
with cTj given by flTTI) in terms of aj. Doing the integral in fl5S]) . we obtain 

^^^^^ =^-^{{l-t ^'a) A^'^^^^) + t 

+ ^ zaAi(^,A)Ai'(2,A)}. {D = 3) (A.23) 

In order to get the explicit expressions for r(xj) or ri(xj), one may apply using 

^P"{x.d = ^Ai\z,x). P = 3) (A.24) 

8m 327r 

Using the expansions (lA.lOp of the Airy function and (]A.4p . we find the leading- 
order oscillating terms in 3D: 

5T{xi) = (A - ajXi) 5p{xi) , (A. 25) 

fulfilling the LVT d^l]), and 

= ^ cos(20a) , (A.26) 

which is by one order h higher than the quantities in flA.25p . 

The densities for D = 5, 7, . . . may be obtained similarly by successive partial 
integrations, but we refrain here from working out the analytical results. Unfortunately, 
we found no simple analytic forms of the densities for even values of D. 
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Appendix B. Explicit densities and relations for the one-dimensional box 



Here we give some explicit results for the one-dimensional box defined in ( H6l) . The 
normalised wave functions fulfilling the Dirichlet boundary condition are 

(j)n{x) = \^2/Lsm{mTx/L) , n= 1,2,3,... (B.l) 

and the eigenvalues are 



E„ 



(B.2) 



2mL2 

The density for particles filling M = N/2 levels (with spin factor 2) becomes (cf. also 

mi 1131111]) 

4 ,^ 

p{x) 



4 2/ If,. sin[(2M + l)7rx/Lll 
- y sin' nvrx/L = - <^ 2M + 1 ^— ,;, ' ^ } 

n=l ^ \ I / y 



2M 

~17 



6p{x) 



(B.3) 



The constant term in the last line is the TF density ptf = 2M/L = N/ L, which can be 
expressed in terms of the Fermi energy Atf by 



Ptf 



N 

T 



2 

TT 



2mA 



TF 



1/2 



A 



TF 



En 



N 



(B.4) 



in agreement with ( 1T0|) for D = 1 and V{v) = 0. The oscillating term in (IB.3P can be 
written as 

Sp{x) = - [2 sin2(M7rx/L) - sin(2M7rx/L) cot(7rx/L)] . (B.5) 

Differentiating this function twice with respect to x, we see that it fulfills, to leading 
order in M, the asymptotic relation 



^ Vas(a;) = 4AtfW 



(B.6) 



This is the equivalent of fl27j) valid asymptotically for IHOs. 
The kinetic-energy density r(x) becomes 

M 



iE, 



TiX 



L 

2Eo 



- sin' {mix / L) 



n=l 
M 



y^ n' [1 — cos(2n7rx/L)] . 



(B.7) 



n=l 



Summing analytically and rearranging terms, we obtain 



2En 



TiX) 



My 3 - M' - sin(2M7rx/L) cot(7ra;/L) 



-sin'(M7rx/L) + 0(M) 
The constant term in the first line is again the TF part: 
2En 



ttf 



L 



2 2m 3/2 



(B.8) 



(B.9) 
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in agreement with (fTTj) for Z) = 1. The leading-order oscillating term in (IB.Sp is 



6tM = ^M' 



sin(2M7rx/L) cot(7rx/L) 



+ sin^(M7rx/L) . (B.IO) 

Combining this with f lB.Sp . it is easy to see that the differential form of the LVT 
derived for IHOs is satisfied here, too, with the proviso V{x) = inside the box: 

St^{x) = \tf5p{x). (B.ll) 

The kinetic-energy density ti{x) becomes 

AE ^ 

ri(x) = J]n2cos2(n7rx/L) . (B.12) 

n=l 

To calculate ^(x), we take the average of (IB. 70 and flB.12p . The sums of squares of sine 
and cosine terms under the summation over n combine to a constant density depending 
only on M, whose asymptotically leading part is the TF kinetic-energy density: 

^ = ^ 1 M{M + 1)(2M + 1) = TTF + C(M2) . (B.13) 

Consequently, the oscillating parts of the two kinetic-energy densities fulfill the relation 
( 123|) obtained for IHOs, replacing the variable r by x: 

6ti{x) = -6t{x) . (B.14) 

The TF functional f lT3|l for the kinetic-energy density for D = 1 is 

ttfIptf] = Ptf ■ (B-15) 

If we insert p(x) from (]B.3P into this functional and expand up to first order in 6p{x), 
we find that the oscillating term is identical with 6tqs{x) given in (IB.lOp . Thus, the TF 
functional relation f lB.lSP holds also for the exact densities of the one-dimensional box 
including the leading-order oscillating terms: 

TTF[pix)] = Ttf + 6tUx) + 0(M) ^ r(x) , (B.16) 

as it was shown in (129|) for IHOs in arbitrary dimensions. 

We should emphasise that, as in the previous examples, the relations (IB. lip and 
(IB.16P do not hold close to the turning points x = and x = L. 

We note that the density oscillations caused by Dirichlet or Neumann boundary 
conditions in one dimension have been interpreted as the manifestation of a "fermionic 
Casimir effect" in [H] (and further references quoted therein). 
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Appendix C. (Integro-) differential equations for the density p{r) 

In this appendix, we briefly discuss some (integro-) differential equations for the density 
p(r) of a system with radial symmetry which are exactly valid for IHOs and linear 
potentials. 

Substituting ( 12T]) into ( IT7|) . we obtain an integro-differential equation (IDE) for the 
spatial density p(r) alone: 

Ap(r) + V{r)p{r) + / V'{q)p{q) dg = Am p(r) . (C. 1) 



8m ' ' ' V /' V / 2 

This is a Schrodinger-type equation, including a non-local potential, with eigenvalue Am 
(Fermi energy). It is exact for IHOs with M flUed shells, using the Fermi energy Am in 
(fT9l) . as shown in [H]. Since the relations (|2T1) and (flTll have been shown in section [331 
to hold also for the liner potential (150]) . the IDE (IC.ip is exact also for this potential, 
provided that r is replaced by any of the Cartesian coordinates Xi. 

Differentiating both sides of ( IC.ip . we can rewrite it as a third-order differential 
equation (30DE) for p{r): 

^ Aap(^) + [A,, - V{r)] -ip(r) + ^V\r)p{r) = . (C.2) 
8m dr dr 2 



This equation had been previously derived for IHOs with D = 1 in |15] and with D = 2 
in [6]. Its form for D = 3 was surmised and numerically tested in [7], and general 
solutions for p{r) in the three-dimensional case were discussed in [TO] . 

For D = 1 dimensional systems, we can expect the IDE (IC.lll to be approximately 
valid, since the SLVT flHT]) is exact and the generalized LVT fl78]) numerically found to 
be well fulflUed everywhere. Therefore, we propose the approximate generalized IDE 
for any differentiable potential V{x): 

- ^ p"(x) + V{x)p{x) + -J V'ix')pix') dx' ^ A pix) , (C.3) 
and the corresponding 30DE: 

^ p"'{x) + [A - V{x)] p'ix) + ^ V'{x)pix) ^ . (C.4) 
8m 2 

The generalization of ( IC.ll) and ( ]C.2I) in D > 1 dimensions poses, however, a 

problem. In the interior region, where (179]) and ( 183]) have to be used without the 

correction terms in brackets {...}, the elimination of ^{r) no longer leads to (integro-) 

differential equations for the density p{r) alone. Taking careful account of the roles of 

the regular and irregular oscillating parts of the density, we would e.g. have to propose 

the following approximate generalized IDE: 

- ^^P{r) + V{r)p{r) + j V'{q)p{q) dq 

^ A[pETF(r) + MO] |+AWW|- (C.5) 
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If the surface correction is included, the full density p(r) appears on the r.h.s. and hence 
the IDE makes sense. In the interior, however, the irregular oscillations Si^-j-lr) are absent 
and we have no longer an IDE for one single function. 




r 



Figure CI. Test of the integro-difFerential equation (jC.5[) without surface correction 
for the three-dimensional potential V{r) = r^/A with N = 42094 (units: h^m=l). 
Solid line: l.h.s., dotted line: r.h.s. of (jC.Sp . 




r 



Figure C2. Same as in figure [CT] but including the surface correction. 



We test (IC.Sp numerically for = 91330 particles in the D = ?> dimensional 
potential V{r) = r^A by comparing both sides with each other. In figure [CTI the surface 
correction is left out. While it fails therefore to reproduce the exponential tail in the 
surface, the equation ( 1C.5I) is seen to very well fulfilled in the interior region. In figure 
C2t the surface correction is included. The quantum-mechanical tail of the density is 
now exactly reproduced, while the error in the interior, which is proportional to Si^plr), 
is still reasonably small. 
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However, as stated above, the equation (IC.Sp without surface correction cannot 
be used to find the full density p(r) for a given smooth potential, since the regular 
oscillating part 6rp{r) is a priori now known. 
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